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The  development  of  a  new  hybrid  CFD  approach  using  fully  coupled  Reynolds-Averaged 
Navier-Stokes  (RANS)  and  Particle-based  Vorticity  Transport  Method  (PVTM)  solvers  is 
summarized  in  this  paper.  The  methodology  is  applied  to  predict  the  performance  and  wake 
parameters  of  an  isolated  wing  in  low  speed  flow  condition.  The  approach  divides  the  flow 
field  into  several  regions  and  uses  appropriate  flow  solvers  according  to  the  dominant 
physical  features  of  the  flow  in  each  region.  The  near  body  flow  field,  which  is  dominated  by 
the  effect  of  viscosity  and  geometry,  is  resolved  using  a  3D  compressible  RANS  solver.  The 
flow  field  outside  of  the  RANS  regions,  which  is  primarily  dominated  by  the  vortices  being 
shed  from  the  aerodynamic  surfaces,  is  simulated  using  a  viscous  implementation  of  the 
Particle-based  Vorticity  Transport  Method.  The  coupling  methodology  and  the  appropriate 
information  being  transferred  between  the  two  solvers  are  also  outlined.  A  stability 
enhancement  procedure  is  implemented  for  the  PVTM  analysis.  The  results  obtained  using 
the  coupled  RANS/PVTM  analysis  compare  well  with  experimental  data,  in  particular  the 
pressure  distribution,  sectional  load,  and  tip  vortex  parameters  (swirl  velocity,  core  location, 
and  core  size). 
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Nomenclature 

strength  of  vortex  particle  ( j  codV) 

vortex  core  function 

density 

kinematic  viscosity 
vorticity  vector 

differential  area  of  outer  boundary 

wing  chord 

lift  force  coefficient 

drag  force  coefficient 

time  step 

differential  Volume 

differential  distance  from  vortex  particle 

vortex  splitting  coefficients 

outward  normal  vector 

pressure 

velocity  vector 

Reynolds  number 

wing  span 

vorticity  Source 

time 

location  vector 
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I.  Introduction 


Accurately  predicting  the  long-term  dynamics  of  the  wake  produced  by  a  wing  under  high  Reynolds  number 
flow  conditions  is  still  one  of  the  most  challenging  tasks  for  CFD  simulations.  In  the  simplest  form,  the  wake 
behind  a  wing  in  these  flow  conditions  usually  consists  of  vortex  sheet  from  the  boundary  layer  and  a  strong  tip 
vortex  from  3D  finite  wing  effects.  The  tip  vortex  is  quickly  formed  and  is  normally  well  organized  before  leaving 
the  wing  trailing  edge1.  A  typical  velocity  profile  of  the  tip  vortex  can  be  found  in  Ref.  2.  The  tip  vortex  continues 
strengthening  after  leaving  the  trailing  edge  by  rolling  in  the  vortex  sheet.  Then  the  viscous  diffusion  effect  causes 
the  tip  vortex  to  grow  in  size  and  lose  strength  very  slowly  as  it  convects  downstream.  The  wake  behind  a  real 
aircraft  is  much  more  complex3.  Normally,  tip  vortices  can  stay  coherent  for  some  distance  behind  the  wing,  and  the 
lifespan  of  the  tip  vortices  is  highly  dependent  on  atmospheric  conditions4.  After  that  the  tip  vortices  can  start 
interacting  with  other  vortices  and  cause  wake  instability5  that  has  been  observed  in  flight  tests6.  This  instability  is 
very  complex  phenomenon7,  and  involves  the  merging  of  unequal  strength  counter-rotating  vortex  pairs8.  The 
simulations  of  such  complex  flow  conditions,  however,  are  beyond  the  scope  of  this  paper. 

Many  researchers  have  successfully  been  simulating  such  wake  flow  using  increasingly  complicated  approaches. 
The  earliest  attempts  simplified  the  wake  structure  and  modeled  it  using  a  vortex  filament  method9  or  vortex  panel 
method10.  A  classical  wing  theory  has  also  been  modified  to  simulate  the  wake  behind  a  wing11.  Recently,  various 
Computational  Fluid  Dynamics  (CFD)  methods  have  been  used  commonly  to  calculate  such  flow,  and  in  one  such 
study,  two  CFD  solvers  (RANS  and  hybrid  Euler/LES  solvers)  were  used  to  compute  aircraft  wakes  very  far 
downstream,  with  an  equivalent  distance  of  2.89  nautical  miles  for  a  full  scale  Airbus  A34012.  In  addition,  vortex 
particle  methods  have  been  used  to  simulate  high  resolution  wake  from  aircraft  using  billions  of  particles13. 

Conventional  CFD  calculations  can  simulate  the  generation  and  evolution  of  the  vortex  sheet  and  tip  vortex  well 
if  sufficient  numbers  of  grid  points  are  used  to  resolve  those  flow  features.  These  calculations  require  localized 
refinement  near  high  velocity  gradient  regions  such  as  the  leading  and  trailing  edges,  boundary  layer,  vortex  sheet, 
and  tip  vortices).  High-order  CFD  algorithms  can  also  be  used  to  improve  the  accuracy  of  the  calculation  in  the  high 
gradient  regions.  Even  with  the  grid  refinement  and  high-order  CFD  algorithms,  the  conventional  CFD  simulations 
still  suffer  from  excessive  numerical  dissipation,  which  typically  causes  the  tip  vortex  to  diffuse  prematurely.  This 
numerical  diffusion  often  causes  the  tip  vortex  core  size  to  grow  faster  than  it  should  based  on  physics.  Though  the 
core  size  becomes  large,  the  vortex  strength  is  preserved,  which  may  be  adequate  for  many  applications.  However, 
for  certain  specific  applications  like  noise  and  vibration  calculations  or  wake  interactions,  accurate  predictions  of  the 
vortex  core  size  and  velocity  profile  are  very  important. 

Recently,  a  new  hybrid  approach  using  fully  coupled  Reynolds-Averaged  Navier-Stokes  (RANS)  and  Particle- 
based  Vorticity  Transport  Method  (PVTM)  solvers  was  introduced  to  simulate  rotorcraft  wake  flow14.  The  approach 
divides  the  flow  field  into  several  regions  and  uses  appropriate  flow  solvers  according  to  the  dominant  physical 
features  of  the  flow  in  each  region.  The  first  region  covers  the  flow  field  near  aerodynamic  surfaces  (extending 
about  1  chord  from  the  surfaces),  where  the  flow  features  are  dominated  by  the  effects  of  boundary  layer  viscosity 
and  the  geometry  of  the  airfoil  and  blade  planform.  The  near  body  flow  field  is  resolved  using  a  3D  compressible 
RANS  solver15.  Outside  of  the  RANS  regions,  the  flow  field  is  primarily  dominated  by  the  vortices  being  shed  from 
the  aerodynamic  surfaces.  This  vortex-dominated  flow  region  is  simulated  using  a  Particle-based  Vorticity 
Transport  Method  (PVTM).  The  influence  of  this  far  field  flow  region  is  transferred  to  the  RANS  regions  using  the 
field  velocity  approach15.  By  modeling  the  far  field  with  PVTM,  the  shed  vorticity  can  remain  well  organized  and  in 
particular  the  tip  vortices  can  maintain  their  compact  and  stable  cores  for  a  longer  period  of  time  than  vortices 
simulated  using  conventional  RANS/CFD.  This  approach  is  particularly  suitable  for  rotorcraft  applications  where 
maintaining  the  correct  velocity  gradient  in  the  vortex  core  is  critical  for  acoustic  and  vibratory  loads  calculations, 
which  depend  greatly  on  the  location,  size,  and  strength  of  the  vortex.  A  previous  study  showed  that  the 
methodology  is  viable,  but  did  not  illustrate  quantitatively  how  well  the  approach  predicted  the  near  body  flow  field 
and  associated  vortex-dominated  flow  field14. 


II.  Focus  of  the  Present  study 

In  the  present  study,  the  coupled  RANS/PVTM  method  is  validated  against  wind  tunnel  test  data.  Comparisons 
with  measured  pressure  distribution,  loadings,  and  vortex  parameters,  and  the  corresponding  results  from  the  full 
RANS  and  coupled  RANS/PVTM  simulations  are  presented  for  a  semi-span  NACA  0015  wing.  The  objectives  are 
to  verify  that  the  coupled  approach  accurately  models  the  wake  and  properly  maintains  the  tip  vortex  geometry. 
Although  RANS  calculations  are  provided  for  reference,  the  pure  RANS  calculations  do  not  take  advantage  of 
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higher  order  or  refinement  techniques,  so  they  do  not  represent  the  state  of  the  art  in  wake  modeling  with  RANS 
CFD.  Rather,  they  represent  a  baseline  to  gauge  the  effectiveness  of  an  alternative  method  to  grid  refinement. 

III.  Approach  and  Method 

The  conventional  RANS  and  coupled  RANS/PVTM  analyses  are  used  to  calculate  the  wing  performance  and 
wake  parameters  for  comparison  with  experimental  data.  The  results  from  the  first  approach  which  simulates  the 
entire  flow  field  using  only  RANS  calculation  are  used  as  a  baseline.  This  baseline  is  compared  to  the  results 
obtained  using  the  hybrid  RANS/PVTM  approach  which  simulates  the  flow  in  the  near  and  far  fields  separately 
using  RANS  and  PVTM,  respectively.  The  following  sections  describe  computational  setup  and  summarize  the 
development  of  the  coupled  RANS/PVTM  methodology. 

A.  Modeling  Experimental  Test  Conditions 

A  comprehensive  experiment  to  measure  the  detailed  progression  of  the  tip  vortex  from  a  semi-span  NACA 
0015  wing16,  conducted  at  the  NASA  Ames  7-  by  10-Foot  Subsonic  Wind  Tunnel,  is  used  as  the  reference  data  to 
validate  the  results  from  the  RANS/PVTM  simulations.  The  experiment  used  a  pressure-instrumented,  untwisted 
semi-span  NACA  0015  wing  with  a  square  tip  (meaning  the  end  is  flat  and  perpendicular  to  the  span  axis).  A  round 
end  cap  was  installed  for  some  conditions  to  change  the  tip  geometry  of  the  wing  from  square  tip  to  rounded  tip. 
The  experiments  were  conducted  for  a  small  range  of  Reynolds  numbers  between  lxlO6  and  3xl06.  Velocity 
profiles  across  the  tip  vortex  were  measured  at  various  downstream  locations  up  to  6  chords  behind  the  wing  using  a 
two-component  laser  velocimeter.  The  data  provided  by  this  test  include  (i)  the  chordwise  pressure  distribution 
along  several  span  locations;  (ii)  the  velocity,  location,  and  core  size  of  the  tip  vortex;  (iii)  integrated  sectional  lift, 
drag,  and  moment  coefficients  at  various  span  locations. 

The  experimental  configuration  and  condition  described  in  Ref.  16  is  used  to  develop  computational  models  of 
the  semi-span  NACA  0015  wing.  The  rectangular  wing  has  a  constant  and  untwisted  NACA  0015  profile  along  the 
span  with  a  full-span  aspect  ratio  of  6.6.  The  wing  was  mounted  on  a  vertical  supporting  end  plate  that  fitted  1  foot 
away  from  the  side  wall  of  the  tunnel.  This  end  plate  prevents  the  formation  of  the  trailing  vortex  near  the 
mounting  point  and  effectively  creates  an  infinite  wing  (2D)  boundary  condition  at  the  mounting  point.  Several 
angles  of  attack  were  tested  in  the  experiment.  For  the  current  study,  only  one  angle  of  attack  of  12°  is  simulated. 
An  additional  0.51°  is  added  to  the  angle  of  attack  in  the  simulations  as  a  correction  for  the  closed-tunnel  wall 
effect16.  The  free-stream  velocity  corresponds  to  a  Mach  number  of  0.13  and  Reynolds  number  of  1.5xl06. 

B.  Full  Reynolds-Averaged  Navier-Stokes  (RANS)  Analysis 

The  baseline  or  full  RANS  simulation  results  in  the  present  study  are  obtained  using  the  TURNS  3D 
compressible  solver15.  The  analysis  solves  the  RANS  equations  for  primitive  variables  p ,  pV ,  pE  inside  a  3D 
structured  grid  using  a  finite  volume  approach.  The  grid  used  in  the  calculation  for  the  square  tip  wing  is  presented 
in  Figure  1  with  the  dimension  of  399x131x97  grid  points.  The  grid  extends  to  about  6  chords  in  all  directions  from 
the  wing  surface,  and  grid  points  are  distributed  a  priori  and  mostly  concentrated  in  regions  with  high  velocity 
gradient  such  as  near  the  leading  edge,  boundary  layer,  tip  vortex,  and  vortex  sheet.  The  boundary  condition  at  the 
root  of  the  wing  is  set  to  be  extrapolated  from  flow  parameters  inside  the  RANS  domain,  simulating  the  infinite 
wing  (2D)  boundary  condition  to  match  the  experiment. 

C.  Coupled  RANS/Particle-based  Vorticity  Transport  Method 

In  the  coupled  RANS/PVTM  analysis14,  the  flow  field  is  divided  into  a  near-body  grid  and  a  far  field  region. 
The  flow  field  in  the  near-body  grid  is  solved  using  the  same  3D  compressible  RANS  solver  and  boundary 
conditions  described  in  the  previous  section,  but  the  domain  is  much  smaller.  The  RANS  domain  extends  to  about  1 
chord  length  in  all  directions  except  behind  the  trailing  edge,  where  it  extends  only  to  about  0.5  chord  length,  see 
Figure  2.  The  small  extent  behind  the  trailing  edge  minimizes  the  dissipation  of  vorticity  to  be  transferred  into 
PVTM  domain.  This  small  RANS  grid  has  a  dimension  of  249x131x79  grid  points.  The  same  2D  boundary 
condition  at  the  wing  root  is  used.  In  addition,  the  induced  velocity  from  the  vortex  particles  in  the  far  field  is 
included  in  the  RANS  calculation  as  field  velocity15. 

Outside  the  near  body  RANS  domain,  the  flow  field  is  represented  by  collections  of  three-dimensional  vortex 
particles  similar  to  those  presented  in  Ref.  17.  Each  vortex  particle  has  two  vector  quantities,  namely  location  and 
strength,  associated  with  it.  The  strength  of  the  vortex  is  a  volume  integration  of  the  vorticity  field  around  the 
particle,  and  the  evolution  of  the  strength  is  governed  by  the  vorticity  transport  equations: 
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(1) 


\^dVt  =  j \co-  Vu]dVi  +  J[vV2«}/\/  +  5. 


Figure  1.  Grid  for  full  RANS  calculation  of  the  square  tip  wing  (5M 
grid  points). 


where  u  is  the  local  velocity,  cois 
vorticity,  v  is  kinematic  viscosity, 
and  Si  is  the  vorticity  source  inside 
a  differential  volume  dVj.  The 
derivation  of  the  vorticity  source 
will  be  explained  in  a  separate 
section  below.  The  location  of  the 
particle  is  changed  due  to  local 
velocity  induced  by  other  vortex 
particles  in  the  flow  field.  The 
detailed  derivation  and 
implementation  of  the  PVTM 
analysis  can  be  found  in  Ref.  14, 
and  is  briefly  summarized  below. 

For  convenience  in  the  PVTM 
domain,  a  multi-level  background 
Cartesian  coordinate  system  shown 
in  Figure  3  is  used  to  facilitate 
fundamental  vortex  particle 
operations  such  as  volume 
integration,  tracking,  merging, 
splitting,  and  induced  velocity 
calculation.  The  PVTM  coordinate 
system  is  different  than  a  grid 
system  in  that  no  variables  are 
being  solved  at  grid  points.  Its 
purpose  is  for  numerical 
integration  and  for  tracking  of 
vortex  particles.  The  coordinate 
system  has  fine  (0.05c)  and  coarse 
(0.2c)  levels  which  extend  to  8c 
and  20c  behind  the  wing, 
respectively.  To  reduce 

computational  cost,  only  one 
particle  is  allowed  in  each  cell.  If 
two  or  more  particles  move  into 
the  same  cell  they  are  merged  with 
the  new  strength  being  the  sum  of 
the  strength  of  the  merged 
particles,  and  the  new  location 
being  the  strength-weighted 
centroid  of  the  particles.  The 
induced  velocity  is  calculated 
using  the  Biot-Savart  law,  and  the 
vortex  stretching  term,  tf>Vw,  is 
calculated  by  differentiating  the  velocity  field.  The  multi-level  particle  mesh  approach  is  employed  to  reduce 
computational  effort  for  calculating  the  induced  velocity  and  stretching  term. 

The  vorticity  source  for  the  PVTM  domain  is  derived  from  the  convection  of  the  vorticity  field  from  the 
boundary  of  the  RANS  domain.  For  each  RANS  cell  that  is  on  the  boundary,  the  vorticity  convected  from  the 
RANS  domain  to  the  PVTM  domain  is  defined  as  follows: 


Figure  2.  Grid  for  coupled  RANS/PVTM  calculation  of  the  square  tip 
wing  (2.5M  grid  points). 
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Figure  3.  A  2D  view  of  a  representative  3D  multi-level 
background  grid  for  PVTM  domain. 


where  ar  is  the  strength  of  the  released 
particle,  Xr  is  the  location  of  the  released 
particle,  A  is  area  of  the  outer  boundary,  h  is 
outward  normal  vector  and  Xmid  is  the  mid¬ 
point  of  the  area  A. 

The  viscous  diffusion  effect  is  simulated 
using  the  vorticity  redistribution  method18,19, 
which  is  performed  separately  from  the 
particle  convection  step.  In  this  study,  a 
single  vortex  particle  is  split  into  five  smaller 
particles  in  the  diffusion  step.  The  locations 
of  the  five  particles  are  defined  by  the 
vertices  of  a  randomly  oriented  tetrahedron 
and  its  center  (see  Figure  4).  The  number 
and  locations  of  the  diffused  particles  are 
somewhat  arbitrary,  and  they  are  chosen  this 
way  to  maintain  the  symmetry  of  the 
vorticity  field.  The  strengths  of  the  particles 
are  determined  by  a  series  expansion  that 
preserves  vorticity  and  its  moments  (i.e.  \d). 


\oxdx ,  and  \ctK2dx).  The  diffusion  equations  are  summarized  in  Eq.  (3)  in  non-dimensionalized  form. 


do) 


(3) 


where 

thatl4 


xv=(Xj  -x,)/h, 

=  1XfAj  =E4^.j 


Re:  Reynolds  number,  f]j  are  determined  by  series  expansion 
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77  =\(OdVh  ^is  the  vortex  core  function,  and  A td  =  NAt.  Due  to  the  high 
Reynolds  number  flow  condition,  the  diffusion  step  is  performed  at  every 
N  time  steps  to  increase  the  effective  diffusion  distance,  h,  and  avoid 
canceling  this  diffusion  effect  with  the  merging  of  the  particles. 

The  process  in  which  each  vortex  particle,  —  \codVt  ,  is  created, 

changed,  and  evolved  with  time  is  simulated  by  a  two-step  procedure.  The 
first  step  is  a  convection  step  which  is  performed  at  every  time  step  and  is 
summarized  in  Eqs.  (4).  Equations  (4)  are  solved  using  4th  order  Runge- 
Kutta  time  integration  algorithm  with  an  adaptive  time  stepping  procedure, 
which  appropriately  adjusts  the  time  step  size  ( dt )  based  on  local  velocity 
and  stress  fields  to  maintain  desired  computational  accuracy. 


Figure  4.  Schematic  of  vorticity 
redistribution  scheme. 
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The  source  term,  an  is  calculated  from  Eq.  (2).  The  superscripts  ( k ,  k-1)  are  the  time  step  indices,  and  the 
subscripts  (/)  are  particle  number  indices.  The  second  diffusion  step  is  performed  at  every  N  time  steps  and  the 
governing  equation  is  given  in  Eq.  (5). 


j 


Figure  5  Coupling  procedure  for  fully  coupled 
RANS/PVTM  analysis. 


A  flowchart  describing  the  coupling  procedure 
between  RANS  and  PVTM  is  shown  in  Figure  5.  At  the 
start  of  any  time  step,  the  RANS  calculation  is  performed 
to  obtain  flow  parameters  that  include  the  induced 
velocity  from  far-field  vortex  particles  as  field  velocity. 
After  the  RANS  calculation  is  completed,  the  velocity 
field  in  the  RANS  domain  is  transformed  into  vorticity 
field,  and  vortex  particles  are  released  into  the  PVTM 
domain  using  Eq.  (2).  Then  these  newly  released 
particles  are  added  to  the  existing  particles  in  the  far  field. 
Each  particle  evolves  and  convects  in  the  far  field,  then  is 
split  (every  N  time  steps)  to  simulate  the  diffusion  effect. 
This  completes  the  calculation  in  the  current  time  step, 
and  the  whole  process  repeats  for  the  next  time  step. 


D.  Tip  Vortex  Core  Parameter  Identification 

The  identification  of  the  tip  vortex  core  parameters, 
namely  size  and  location,  is  adopted  from  the  procedure 
presented  in  Ref.  16.  In  the  experiment,  an  approximate 
location  of  the  tip  vortex  at  some  distance  behind  the 
trailing  edge  was  measured  using  a  vortex  meter.  Then  a 
laser  velocimeter  was  used  to  measure  swirl  velocity 
across  the  vortex  core  in  the  spanwise 
direction.  The  vortex  core  size  reported  is  the 
distance  between  the  locations  of  the 
maximum  and  minimum  swirl  velocity  across 
the  core.  The  precise  location  of  the  core  was 
determined  to  be  the  mid-point  between  the 
minimum  and  maximum  swirl  velocity. 

Similar  methodologies  are  employed  for 
the  computational  results.  To  identify  the 
vortex  core  parameters  inside  a  RANS 
domain,  the  approximate  location  of  the  tip 
vortex  core  is  determined  by  searching  for  the 
maximum  amplitude  of  vorticity  in  a  given 
plane  behind  the  trailing  edge.  Once  this  is 
known,  the  swirl  velocity  across  the  vortex 
core  is  calculated  by  linear  interpolation,  and 
then  the  core  size  and  precise  location  of  the 
core  are  calculated  in  the  same  manner  as  in 
the  experiment. 

A  slightly  different  procedure  is  used  for 


Figure  6.  Vorticity  profile  calculated  using  coupled 
RANS/PVTM  method  showing  velocity  planes  and  vortex 
particles  in  trailed  wake. 
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identifying  the  vortex  core  parameters  for  the  PVTM  domain,  since  the  PVTM  analysis  does  not  calculate  the 
velocity  field  directly.  At  the  desired  distance  behind  the  trailing  edge,  a  properly  located  velocity  measurement 
plane  is  created  to  record  the  induced  velocity  from  the  vortex  particle  field.  These  are  shown  in  Figure  6  for  1,  2,  4, 
and  6  chords  behind  the  airfoil.  The  induced  velocity  field  from  the  vortex  particles  is  calculated  in  each 
measurement  plane,  similar  to  particle  image  velocimetry  measurements  in  an  experiment.  The  velocity  fields  are 
processed  into  vorticity  with  the  strength  vector  perpendicular  to  the  measurement  plane.  The  vorticity  field  is 
searched  to  find  the  location  of  the  maximum  vorticity,  to  obtain  the  location  of  the  vortex  core.  The  swirl  velocity 
can  be  extracted  directly  from  the  measurement  plane,  to  obtain  the  core  size  in  the  same  way  as  for  the 
measurement  and  RANS -calculated  velocity  field.  For  the  PVTM  vortex  parameters,  200  time  steps  are  averaged  to 
obtain  the  vortex  parameters.  This  is  necessary  because  the  induced  velocities  from  the  particles  are  highly 
dependent  on  the  proximity  of  the  particle  to  the  point  where  the  velocity  is  being  calculated.  Averaging  smoothes 
out  this  nonlinearity  to  obtain  a  more  consistent  result. 

E.  Stability  Enhancement  Procedure  for  Particle  Vortex  Method 

One  of  the  stability  conditions  cited  in  Ref.  17  was  that  the  flow  field  must  be  divergence  free  at  all  time.  This 
requires  the  vorticity  field  to  be  uniformly  distributed  and  align  itself  properly.  Several  modifications  to  the  induced 
velocity  and  stretching  term  calculations  are  made  to  achieve  this  goal  as  shown  below.  The  induced  velocity 
equation  is  changed  from  Eq.  (6a)  to  (6b).  An  attraction  term  (2nd  term)  is  added  to  counteract  the  tendency  of  the 
particles  to  slowly  drift  apart  in  the  simulation,  which  helps  the  vortex  core  to  maintain  its  compactness  and 
coherence. 
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The  stretching  term  is  modified  from  Eq.  (7a)  to  (7b)  to  help  align  the  strength  vector  of  each  vortex  particle 
based  on  the  strength  of  the  surrounding  particle  field.  In  effect,  this  is  analogous  to  the  Particle  Strength  Exchange 
(PSE)  method20. 
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The  coefficients  r|1  and  r\2  control  the  overall  stability  of  the  PVTM  calculation.  They  depend  on  a  number  of 
parameters  including  time  step  size  ( dt ),  and  in  this  particular  case  they  are  chosen  to  be  0.005  and  5,  respectively. 

IV.  Results 


A.  Pressure  Distribution 

Obtaining  the  pressure  distribution  around  the  airfoil  from  the  RANS  calculation  is  straightforward.  The  test 
data  includes  pressure  distributions  at  several  span  stations  along  the  wing.  The  first  thing  to  be  verified  is  the 
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Aerodynamic  Coefficient 


(a)  0.974R 


(b)  0.944R 


(c)  0.845R 


(d)  0.773R 


Figure  7.  Airfoil  surface  pressure  distributions  near  the  wing  tip. 


Figure  8.  Lift  and  drag  coefficients  along  wing  span 
derived  from  surface  pressures. 
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prediction  of  the  pressure  field  around 
the  airfoil  which  is  highly  affected  by  the 
near  field  wake  as  well  as  the  far  field 
wake.  For  the  coupled  method,  the 
influence  of  the  trailed  vortex  must  be 
correctly  accounted  for  inside  the  RANS 
domain  for  accurate  prediction  of  the 
pressure  field  on  the  airfoil  surface. 

Surface  pressures  for  four  radial 
stations  near  the  wing  tip  are  shown  in 
Figure  7.  The  surface  pressures  are 
visibly  affected  by  the  tip  vortex 
outboard  of  0.845,  but  inboard  of  0.845, 
the  pressure  distributions  are  the  same. 
The  plots  show  that  the  hybrid 
RANS/PVTM  method  captures  the 
variation  in  surface  pressure  at  least  as 
well  as  the  baseline  RANS  calculation. 
For  the  outermost  radial  station,  the 
result  from  the  hybrid  RANS/PVTM 
calculation  is  closer  to  the  test  data  near 
the  trailing  edge  (0.8c).  The  surface 
pressures  are  integrated  to  obtain  lift  and 
drag  along  the  wing  span,  shown  in 
Figure  8.  The  span  wise  variation  in  lift 
and  drag  due  to  the  tip  vortex  is  well 
captured  for  both  the  RANS  and  hybrid 
RANS/PVTM  calculations.  These  results 
show  that  the  interchange  of  vorticity  and 
velocity  between  the  near-body  RANS 
grid  and  the  off-body  PVTM  domain  is  properly 
implemented. 


B.  Vorticity  Sources 

Using  the  near-body  RANS  grid  can  yield 
several  advantages  over  some  simplified 
methods  such  as  lifting  line  theory  or  vortex 
panel  method.  The  RANS  calculations  provide 
much  more  detailed  flow  field  information  such 
as  pressure,  loading,  and  shape  of  the  tip  vortex, 
than  other  methods.  This  is  shown  with  the 
different  profiles  of  the  vorticity  in  the  region  of 
the  tip  vortex  at  the  interface  between  the 
RANS  and  PVTM  domains  for  the  square  and 
rounded  tip  wings  in  Figure  9.  The  vortex  sheet 
being  trailed  behind  the  wing  and  the  tip  vortex 
are  visible,  and  the  different  character  of  the  tip 
vortex  is  evident  between  the  two  tip 
geometries.  For  the  square  tip  wing,  the  tip 
vortex  is  formed  from  smaller  vortices  being 
shed  from  the  upper  and  lower  edges,  while  for 
the  round  tip  wing,  a  single  compact  tip  vortex 
is  shed  from  the  tip.  This  sort  of  detail  is  not 
captured  unless  a  viscous,  3D  RANS  simulation 
is  used  in  the  near-body  region. 
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Square  Tip 


Round  Tip 


Figure  9.  Vorticity  distribution  in 
vortex  sheet  and  tip  vortex  for  square 
and  round  tip  airfoils. 


Downstream  Location  (Chord  Lengths) 


Figure  10.  Vortex  core  size  downstream  of  the  wing. 


C.  Tip  Vortex  Geometry 

The  next  comparison  verifies  the  vortex  geometry,  namely  size  and  location,  of  the  tip  vortex.  Figure  10  shows 
the  vortex  core  size  as  a  function  of  downstream  location  and  Figure  11  shows  its  location.  Both  the  RANS  and 
hybrid  RANS/PVTM  calculations  show  an  offset  in  core  size  near  the  trailing  edge.  Downstream,  however,  the  core 
size  steadily  increases  in  the  RANS  calculation,  while  the  hybrid  calculation  maintains  a  nearly  constant  core  size 
which  compares  well  with  the  test  data. 


Downstream  Location  (Chord  Lengths)  Downstream  Location  (Chord  Lengths) 


(a)  Spanwise  location  (b)  Vertical  location 

Figure  11.  Spanwise  and  vertical  locations  of  tip  vortex  downstream  of  the  wing. 
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Downstream  Location  (Chord  Lengths) 
Figure  12.  Tip  vortex  geometry  calculated  by 
the  RANS/PVTM  method  with  and  without 
additional  attraction  term. 


Spanwise  Distance  from  Tip  (mm) 


Figure  13.  Swirl  velocity  profile  across  tip  vortex 
four  chords  behind  trailing  edge. 


The  location  of  the  tip  vortex  is  shown  in  Figure  1 1 .  The  variation  in  spanwise  location  is  captured  well  in  both 
the  RANS  and  hybrid  RANS/PVTM  methods.  The  correlation  of  vertical  location  is  not  as  good,  but  the  magnitude 
of  the  variation  is  also  only  half  that  of  the  spanwise  location,  making  the  differences  appear  larger.  The  vortex  does 
gradually  drift  upward  like  the  experiment  for  both  cases.  The  vortex  attraction  term  can  be  used  to  improve  the 
correlation.  The  PVTM  results  with  and  without  the  attraction  term  are  compared  to  the  test  data  in  Figure  12.  With 
the  attraction  term,  the  offset  in  core  size  is  removed  such  that  the  calculated  core  size  matches  the  experimental 
data  very  closely.  The  vertical  location  of  the  tip  vortex  is  also  improved. 


D.  Velocity  Profiles 

Velocity  profiles  at  four  and  six  chords  behind  the  airfoil  are  shown  in  Figure  13  and  Figure  14,  respectively. 
For  both  the  RANS  and  hybrid  RANS/PVTM  calculations,  the  peak  velocities  are  under-predicted.  At  6  chords,  the 


numerical  diffusion  in  the  pure  RANS  calculation  i 


Spanwise  Distance  from  Tip  (mm) 


evident  while  the  hybrid  method  maintains  a  stable  and 
compact  core  size.  The  attraction  term  results  in  a  much 
sharper  gradient  through  the  core  and  hence  a  smaller  tip 
vortex.  From  these  figures,  it  is  not  as  obvious  what  the 
core  size  is  as  in  Figure  12,  for  example,  but  the  general 
trend  can  be  established.  It  is  clear  that  the  attraction 
term  causes  the  peak  velocity  to  be  over-predicted,  so 
careful  selection  of  the  attraction  parameters  is 
important. 

Finally,  the  effects  of  viscous  diffusion  are  shown  in 
Figure  15.  The  three  calculations  are  the  baseline 
PVTM,  PVTM  with  the  attraction  term,  and  PVTM  with 
both  the  attraction  term  and  viscous  diffusion.  It  is  clear 
that  the  diffusion  makes  no  difference  to  the  resolution 
of  the  plot.  This  is  not  entirely  unexpected  as  6  chords  is 
a  very  young  wake  age  and  the  effects  of  real  viscosity  in 
the  tip  vortex  should  not  be  evident  until  much  farther 
downstream. 

V.  Conclusions 


Figure  14.  Swirl  velocity  profile  in  tip  vortex  six  A  hybrid  Reynolds- Averaged  Navier-Stokes  and 

chords  behind  trailing  edge.  Particle-based  Vorticity  Transport  Method  solution 
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Spanwise  Distance  from  Tip  (mm) 
Figure  15.  Velocity  profile  in  tip  vortex  for 
PVTM  calculation  six  chords  behind  the  trailing 
edge,  showing  influence  of  attraction  term  and 
viscous  diffusion. 


methodology  has  been  presented  and  validated  against  test 
data  for  a  fixed  wing  in  a  wind  tunnel.  Based  on  the 
results  observed  in  this  study,  the  following  conclusions 
are  offered. 

The  combination  of  vortex  shedding  from  the  near 
body  domain  into  the  PVTM  domain  and  the  field 
velocity  method  for  communicating  the  influence  of  the 
far  wake  on  the  near  body  domain  is  an  effective  coupling 
methodology.  The  three-dimensional  effects  of  the  wing 
tip  were  properly  modeled  in  the  near  body  domain. 
Correlation  between  calculated  lift  and  drag  coefficients 
and  measured  data  was  satisfactory. 

The  hybrid  method  resulted  in  less  numerical 
dissipation  of  the  tip  vortex  than  the  baseline  RANS 
calculation,  which  showed  significant  diffusion  by  4-6 
chords  wake  age.  For  the  hybrid  method,  the  core  size 
was  initially  larger  than  in  the  experiment,  but  maintained 
nearly  constant  size,  similar  to  the  test  data.  The  location 
of  the  vortex  was  satisfactorily  captured  by  both  the  pure 
RANS  and  the  hybrid  method. 

An  added  attraction  term  improved  correlation  with 
test  data  by  reducing  the  size  of  the  vortex  core.  Viscous 
diffusion  did  not  have  an  effect  on  the  results  for  wake 
ages  up  to  6  chords. 
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